pkgs <- c(
  "foreign", "dplyr", "ggplot2", "stargazer", "stringr",
  "lubridate", "readstata13", "xtable", "haven","lmtest"
)
bools <- pkgs %in% installed.packages()[,1]
to_install <- NA
if(mean(bools) != 1){
  to_install <- pkgs[!bools]
  for(i in 1:length(to_install)){
    install.packages(
      to_install[i],
      dependencies = T)
  }
}

###MODIFY:
###
#load packages, remove objects, load data
sapply(pkgs, require, character.only = T); rm(pkgs); rm(to_install); rm(bools)
data <- read_dta("final_dataset.dta")

####populate w new variables needed
data$authority <- data$authority2
data$illiterate <- ifelse(data$P1_22 ==2 | data$P1_23==2, 1, 0)
data$indigenous <- ifelse(data$P1_20 == 1, 1, 0)

data$sentenced <- NA
data$sentenced[data$P5_3 == 1] <- 1
data$sentenced[data$P5_3 == 2] <- 2
data$sentenced[data$P5_3 == 3] <- 3

data$cat_age <- NA
data$cat_age[data$age_at_arrest <= 25] <- 1
data$cat_age[data$age_at_arrest <= 35&data$age_at_arrest > 25] <- 2
data$cat_age[data$age_at_arrest <= 45&data$age_at_arrest > 35] <- 3
data$cat_age[data$age_at_arrest <= 55&data$age_at_arrest > 45] <- 4
data$cat_age[data$age_at_arrest <= 65&data$age_at_arrest > 55] <- 5
data$cat_age[data$age_at_arrest > 65] <- 6


# APPENDIX TABLE A6 --------------------------------------------------------
###Covariates ebfore and after reform

pre <- c(
  t.test(data$SEXO[data$allreform_with_fed == 0], data$SEXO[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$illiterate[data$allreform_with_fed == 0], data$illiterate[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$indigenous[data$allreform_with_fed == 0], data$indigenous[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$authority[data$allreform_with_fed == 0] == 1, data$authority[data$allreform_with_fed == 1] == 1)$estimate[1], 
  t.test(data$authority[data$allreform_with_fed == 0] == 2, data$authority[data$allreform_with_fed == 1] == 2)$estimate[1], 
  t.test(data$authority[data$allreform_with_fed == 0] == 3, data$authority[data$allreform_with_fed == 1] == 3)$estimate[1], 
  t.test(data$authority[data$allreform_with_fed == 0] == 4, data$authority[data$allreform_with_fed == 1] == 4)$estimate[1], 
  t.test(data$authority[data$allreform_with_fed == 0] == 5, data$authority[data$allreform_with_fed == 1] == 5)$estimate[1], 
  t.test(data$theft[data$allreform_with_fed == 0], data$theft[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$homicide[data$allreform_with_fed == 0], data$homicide[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$kidnap[data$allreform_with_fed == 0], data$kidnap[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$illegal_weapons[data$allreform_with_fed == 0], data$illegal_weapons[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$rape[data$allreform_with_fed == 0], data$rape[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$drug_possess[data$allreform_with_fed == 0], data$drug_possess[data$allreform_with_fed == 1])$estimate[1], 
  t.test(data$drug_commerce[data$allreform_with_fed == 0], data$drug_commerce[data$allreform_with_fed == 1])$estimate[1],
  t.test(data$sentenced[data$allreform_with_fed == 0] == 1, data$sentenced[data$allreform_with_fed == 1] == 1)$estimate[1],
  t.test(data$sentenced[data$allreform_with_fed == 0] == 2, data$sentenced[data$allreform_with_fed == 1] == 2)$estimate[1],
  t.test(data$sentenced[data$allreform_with_fed == 0] == 3, data$sentenced[data$allreform_with_fed == 1] == 3)$estimate[1]
)
post <- c(
  t.test(data$SEXO[data$allreform_with_fed == 0], data$SEXO[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$illiterate[data$allreform_with_fed == 0], data$illiterate[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$indigenous[data$allreform_with_fed == 0], data$indigenous[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$authority[data$allreform_with_fed == 0] == 1, data$authority[data$allreform_with_fed == 1] == 1)$estimate[2], 
  t.test(data$authority[data$allreform_with_fed == 0] == 2, data$authority[data$allreform_with_fed == 1] == 2)$estimate[2], 
  t.test(data$authority[data$allreform_with_fed == 0] == 3, data$authority[data$allreform_with_fed == 1] == 3)$estimate[2], 
  t.test(data$authority[data$allreform_with_fed == 0] == 4, data$authority[data$allreform_with_fed == 1] == 4)$estimate[2], 
  t.test(data$authority[data$allreform_with_fed == 0] == 5, data$authority[data$allreform_with_fed == 1] == 5)$estimate[2], 
  t.test(data$theft[data$allreform_with_fed == 0], data$theft[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$homicide[data$allreform_with_fed == 0], data$homicide[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$kidnap[data$allreform_with_fed == 0], data$kidnap[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$illegal_weapons[data$allreform_with_fed == 0], data$illegal_weapons[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$rape[data$allreform_with_fed == 0], data$rape[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$drug_possess[data$allreform_with_fed == 0], data$drug_possess[data$allreform_with_fed == 1])$estimate[2], 
  t.test(data$drug_commerce[data$allreform_with_fed == 0], data$drug_commerce[data$allreform_with_fed == 1])$estimate[2],
  t.test(data$sentenced[data$allreform_with_fed == 0] == 1, data$sentenced[data$allreform_with_fed == 1] == 1)$estimate[2],
  t.test(data$sentenced[data$allreform_with_fed == 0] == 2, data$sentenced[data$allreform_with_fed == 1] == 2)$estimate[2],
  t.test(data$sentenced[data$allreform_with_fed == 0] == 3, data$sentenced[data$allreform_with_fed == 1] == 3)$estimate[2]
)
p <- c(
  t.test(data$SEXO[data$allreform_with_fed == 0], data$SEXO[data$allreform_with_fed == 1])$p.value, 
  t.test(data$illiterate[data$allreform_with_fed == 0], data$illiterate[data$allreform_with_fed == 1])$p.value, 
  t.test(data$indigenous[data$allreform_with_fed == 0], data$indigenous[data$allreform_with_fed == 1])$p.value, 
  t.test(data$authority[data$allreform_with_fed == 0] == 1, data$authority[data$allreform_with_fed == 1] == 1)$p.value, 
  t.test(data$authority[data$allreform_with_fed == 0] == 2, data$authority[data$allreform_with_fed == 1] == 2)$p.value, 
  t.test(data$authority[data$allreform_with_fed == 0] == 3, data$authority[data$allreform_with_fed == 1] == 3)$p.value, 
  t.test(data$authority[data$allreform_with_fed == 0] == 4, data$authority[data$allreform_with_fed == 1] == 4)$p.value, 
  t.test(data$authority[data$allreform_with_fed == 0] == 5, data$authority[data$allreform_with_fed == 1] == 5)$p.value, 
  t.test(data$theft[data$allreform_with_fed == 0], data$theft[data$allreform_with_fed == 1])$p.value, 
  t.test(data$homicide[data$allreform_with_fed == 0], data$homicide[data$allreform_with_fed == 1])$p.value, 
  t.test(data$kidnap[data$allreform_with_fed == 0], data$kidnap[data$allreform_with_fed == 1])$p.value, 
  t.test(data$illegal_weapons[data$allreform_with_fed == 0], data$illegal_weapons[data$allreform_with_fed == 1])$p.value, 
  t.test(data$rape[data$allreform_with_fed == 0], data$rape[data$allreform_with_fed == 1])$p.value, 
  t.test(data$drug_possess[data$allreform_with_fed == 0], data$drug_possess[data$allreform_with_fed == 1])$p.value, 
  t.test(data$drug_commerce[data$allreform_with_fed == 0], data$drug_commerce[data$allreform_with_fed == 1])$p.value,
  t.test(data$sentenced[data$allreform_with_fed == 0] == 1, data$sentenced[data$allreform_with_fed == 1] == 1)$p.value,
  t.test(data$sentenced[data$allreform_with_fed == 0] == 2, data$sentenced[data$allreform_with_fed == 1] == 2)$p.value,
  t.test(data$sentenced[data$allreform_with_fed == 0] == 3, data$sentenced[data$allreform_with_fed == 1] == 3)$p.value
)

table <- 
  data.frame(
    all_data_pre = pre,
    all_data_post = post,
    all_p = p
  )

xtable(
  table,
  digits = 3
)

# APPENDIX TABLE A12 ------------------------------------------------------
##Table of announcement dates
states <- unique(data$Ent_arrest[data$fecha_publicacion_primera!= ""]) %>% sort
date <- c()
stateout <- c()
rec <- 0
for(i in 1:length(states)){
  tmp <- unique(
    na.omit(data$fecha_publicacion_primera[data$Ent_arrest == states[i]])
  )
  tmp <- tmp[tmp != ""]
  if(length(tmp) > 1){
    print(tmp)
    print(i)
  }
  
    date <- c(
      date, tmp
  )
}
tab <- data.frame(states = states, date = date)
print(xtable(tab), include.rownames=FALSE)

# MAIN TEXT: TABLE 2 -----------------------------------------------------------------

##AVERAGES
before <- c("beat_before", "beat_w_objects_before",
  "crush_before_mpj", "suffocation_before_mpj", "burns_before_mpj", 
  "electricity_before_mpj", "stab_before_mpj",
  "false_charges_threat_before", "harm_family_threat_before",
  "incomunicado_b", "strip_b", "tie_b", "blindfold_b"
  )
pm <- c("pubmin_beat", "pubmin_beat_w_objects",
        "crush_at_pubmin", "suffocation_at_pubmin", "burns_at_pubmin",
        "electricity_at_pubmin", "stab_at_pubmin", 
        "pubmin_false_charge_threat", "pubmin_family_threat",
        "incomunicado_pm", "strip_pm", "tie_pm", "blindfold_pm"
)
rnam <- c("beat", "beat_w_objects",
          "crush", "suffocation", "burns",
          "electricity", "stab", 
          "false_charge_threat", "family_threat",
          "incomunicado", "strip", "tie", "blindfold"
)

rows <- matrix(nrow = 13, ncol = 3)

for(i in 1:13){
  mini_dat <- data[,c(before[i], pm[i])]
  for(j in 1:2){
    mini_dat[is.na(mini_dat[,j]),j] <- 0
  }
  mini_dat$both <- ifelse((mini_dat[,1] + mini_dat[,2]) == 2, 1, 0)
  row <- c(sum(mini_dat[,1], na.rm = T)/nrow(data), 
           sum(mini_dat[,2], na.rm = T)/nrow(data), 
           sum(mini_dat[,3], na.rm = T)/nrow(data))
  rows[i,1:3] <- row
}
rownames(rows) <- rnam
colnames(rows) <- c("Before MP", "At MP", "Both")
rows <- rows*100
xtable(rows, digits = 2)


##COUNTS
rows <- matrix(nrow = 13, ncol = 3)
for(i in 1:13){
  mini_dat <- data[,c(before[i], pm[i])]
  for(j in 1:2){
    mini_dat[is.na(mini_dat[,j]),j] <- 0
  }
  mini_dat$both <- ifelse((mini_dat[,1] + mini_dat[,2]) == 2, 1, 0)
  row <- c(paste0('(', as.character(sum(
    mini_dat[,1], na.rm = T)), ')'), 
           paste0('(', as.character(sum(
             mini_dat[,2], na.rm = T)), ')'), 
           paste0('(', as.character(sum(
             mini_dat[,3], na.rm = T)), ')'))
  rows[i,1:3] <- row
}
rownames(rows) <- rnam
colnames(rows) <- c("Before MP", "At MP", "Both")

xtable(rows)


# TABLE 5 MAIN TEXT -------------------------------------------------------
outcomes <- c("brute_force_dum", "institutional_torture_dum", "threat_dum")

mean_before <- c()
mean_after <- c()
p <- c()

for(i in 1:3){
  mini_dat <- data %>% filter(federal_prison == 0) %>% data.frame()
  mini_dat <- mini_dat[,c("allreform", outcomes[i])]
  mean_before <- c(mean_before, 
                   t.test(mini_dat[mini_dat$allreform == 0, 2],
                          mini_dat[mini_dat$allreform == 1, 2])$estimate[1]
  )
  mean_after <- c(mean_after, 
                  t.test(mini_dat[mini_dat$allreform == 0, 2],
                         mini_dat[mini_dat$allreform == 1, 2])$estimate[2]
  )
  p <- c(p,
         t.test(mini_dat[mini_dat$allreform == 0, 2],
                mini_dat[mini_dat$allreform == 1, 2])$p.value
  )
}
means
mean_after
mean_before

table <- data.frame(
  before_ref = mean_before,
  after_ref = mean_after,
  difference = mean_after - mean_before,
  p = p
)
xtable(table*100,
       digits = 2)



# TABLE 4: Abuse by conviction status -------------------------------------
#nosent: 1 
#sent: 3
#part sent: 2

brf <- c(
  mean(data$brute_force_dum[data$P5_3 == 1], na.rm = T),
  mean(data$brute_force_dum[data$P5_3 == 3], na.rm = T),
  mean(data$brute_force_dum[data$P5_3 == 2], na.rm = T)
)
inst <- c(
  mean(data$institutional_torture_dum[data$P5_3 == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$P5_3 == 3], na.rm = T),
  mean(data$institutional_torture_dum[data$P5_3 == 2], na.rm = T)
)

thr <- c(
  mean(data$threat_dum[data$P5_3 == 1], na.rm = T),
  mean(data$threat_dum[data$P5_3 == 3], na.rm = T),
  mean(data$threat_dum[data$P5_3 == 2], na.rm = T)
)

n <- c(
  sum(data$P5_3 == 1, na.rm = T),
  sum(data$P5_3 == 3, na.rm = T),
  sum(data$P5_3 == 2, na.rm = T)
)

n_abuse <- c(
  sum(data$brute_force_dum, na.rm = T),
  sum(data$institutional_torture_dum, na.rm = T),
  sum(data$threat_dum, na.rm = T),
  NA
)

tab <- cbind(brf, inst, thr, n)
tab <- rbind(tab, n_abuse)

rownames(tab) <- c("Not Sentenced", "Sentenced", "Partly Sentenced", "(Num)")
colnames(tab) <- c("Brute Force", "Inst. Torture", "Threats", "(Num)")
xtable(tab,digits = 4)


# APPENDIX ----------------------------------------------------------------
####Table: implementation dates
#load dates
dat <- read.csv("C:/Users/User/Box/Paper Tortura/datasets/muni merged dates restricted3.csv")
#Extract relevant state codes
#Convert dates to date to allow sorting by time
dat$Implementacion <- dat$Implementacion %>% mdy()
states <- sort(
  unique(
    data$Ent_arrest[!is.na(data$Implementacion)]
  )
)
#Load state names
key <- read.csv("C:/Users/User/Box/Paper Tortura/datasets/untitled folder/ARCH245.csv")
#remove information about municipalities, get rid of duplicates
key <- key %>% dplyr::select(
  CVE_ENT, NOM_ENT
) %>% filter(
  !duplicated(CVE_ENT)
)
#create empty dataframe

dates <- matrix(nrow = length(states), ncol = 2) %>% data.frame
colnames(dates) <- c("state", "dates")
state <- NA
for(i in 1:length(states)){
  #get state code i, record it; record state name
  state <- states[i]
  dates$state[i] <- as.character(key$NOM_ENT[key$CVE_ENT == state])
  #retrieve a sorted vector containing all unique dates for state i,
  #collapse into a single string and store
  dates$dates[i] <- paste(
    sort(
      unique(
        dat$Implementacion[dat$Ent_arrest == state]
      )
    ),
    collapse = ", "
  )
}

print(
  xtable(
    dates
  ), include.rownames=FALSE
)
#How many unique implementation dates are there?
length(unique(dat$Implementacion))



# Table 3: Covariates of torture -----------------------------------------


data$illiterate <- 0
data$illiterate[data$noread_write == 2] <- 1
data$indigenous <- 0
data$indigenous[data$P1_20 == 1] <- 1

data$occupation <- NA
data$occupation[!is.na(data$P2_9)] <- "Other"
data$occupation[data$P2_9 == 4|data$P2_9 == 5] <- "Merchant"
data$occupation[data$P2_9 == 8|data$P2_9 == 6|data$P2_9 == 7|data$P2_9 == 9] <- "Public Security"
data$occupation[data$P2_9 == 11] <- "Private security"
data$occupation[data$P2_9 == 12] <- "Rural worker"
data$occupation[data$P2_9 == 13] <- "Craftsman"
data$occupation[data$P2_9 == 14] <- "Blue collar worker"

data$wealth_q <- NA
data$wealth_q[data$index_wealth < quantile(data$index_wealth, prob = .25, na.rm = T)] <- "0% - 25%"
data$wealth_q[data$index_wealth >= quantile(data$index_wealth, prob = .25, na.rm = T) & data$index_wealth < quantile(data$index_wealth, prob = .5, na.rm = T)] <- "25% - 50%"
data$wealth_q[data$index_wealth >= quantile(data$index_wealth, prob = .5, na.rm = T) & data$index_wealth < quantile(data$index_wealth, prob = .75, na.rm = T)] <- "50% - 75%"
data$wealth_q[data$index_wealth >= quantile(data$index_wealth, prob = .75, na.rm = T)] <- "75%+"

data$cat_age <- NA
data$cat_age[data$age_at_arrest <= 25] <- 1
data$cat_age[data$age_at_arrest <= 35&data$age_at_arrest > 25] <- 2
data$cat_age[data$age_at_arrest <= 45&data$age_at_arrest > 35] <- 3
data$cat_age[data$age_at_arrest <= 55&data$age_at_arrest > 45] <- 4
data$cat_age[data$age_at_arrest <= 65&data$age_at_arrest > 55] <- 5
data$cat_age[data$age_at_arrest > 65] <- 6

data$education <- NA
data$education[data$P1_24_N == 1|data$P1_24_N == 2] <- 1
data$education[data$P1_24_N == 3|data$P1_24_N == 4|data$P1_24_N == 5] <- 2
data$education[data$P1_24_N == 6|data$P1_24_N == 7] <- 3
data$education[data$P1_24_N == 8|data$P1_24_N == 9] <- 4



col1 <- c(
  mean(data$brute_force_dum[data$P3_2 == 1], na.rm = T),
  mean(data$brute_force_dum[data$P3_2 == 2], na.rm = T),
  mean(data$brute_force_dum[data$P3_2 == 4], na.rm = T),
  mean(data$brute_force_dum[data$P3_2 == 3], na.rm = T),
  mean(data$brute_force_dum[data$P3_2 == 5|data$P3_2 == 6], na.rm = T),
  mean(data$brute_force_dum[data$federal_prison == 0], na.rm = T),
  mean(data$brute_force_dum[data$federal_prison == 1], na.rm = T),
  mean(data$brute_force_dum[data$homicide == 1], na.rm = T),
  mean(data$brute_force_dum[data$rape == 1], na.rm = T),
  mean(data$brute_force_dum[data$theft == 1], na.rm = T),
  mean(data$brute_force_dum[data$drug_commerce == 1], na.rm = T),
  mean(data$brute_force_dum[data$drug_possess == 1], na.rm = T),
  mean(data$brute_force_dum[data$extortion == 1], na.rm = T),
  mean(data$brute_force_dum[data$illegal_weapons == 1], na.rm = T),
  mean(data$brute_force_dum[data$kidnap == 1], na.rm = T),
  mean(data$brute_force_dum[data$SEXO == 1], na.rm = T),
  mean(data$brute_force_dum[data$SEXO == 0], na.rm = T),
  mean(data$brute_force_dum[data$education == 1], na.rm = T),
  mean(data$brute_force_dum[data$education == 2], na.rm = T),
  mean(data$brute_force_dum[data$education == 3], na.rm = T),
  mean(data$brute_force_dum[data$education == 4], na.rm = T),
  mean(data$brute_force_dum[data$indigenous == 1], na.rm = T),
  mean(data$brute_force_dum[data$indigenous == 0], na.rm = T),
  mean(data$brute_force_dum[data$cat_age == 1], na.rm = T),
  mean(data$brute_force_dum[data$cat_age == 2], na.rm = T),
  mean(data$brute_force_dum[data$cat_age == 3], na.rm = T),
  mean(data$brute_force_dum[data$cat_age == 4], na.rm = T),
  mean(data$brute_force_dum[data$cat_age == 5], na.rm = T),
  mean(data$brute_force_dum[data$cat_age == 6], na.rm = T),
  mean(data$brute_force_dum, na.rm = T),
  sum(data$brute_force_dum, na.rm = T)
  
)

col2 <- c(
  mean(data$institutional_torture_dum[data$P3_2 == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$P3_2 == 2], na.rm = T),
  mean(data$institutional_torture_dum[data$P3_2 == 4], na.rm = T),
  mean(data$institutional_torture_dum[data$P3_2 == 3], na.rm = T),
  mean(data$institutional_torture_dum[data$P3_2 == 5|data$P3_2 == 6], na.rm = T),
  mean(data$institutional_torture_dum[data$federal_prison == 0], na.rm = T),
  mean(data$institutional_torture_dum[data$federal_prison == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$homicide == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$rape == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$theft == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$drug_commerce == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$drug_possess == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$extortion == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$illegal_weapons == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$kidnap == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$SEXO == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$SEXO == 0], na.rm = T),
  mean(data$institutional_torture_dum[data$education == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$education == 2], na.rm = T),
  mean(data$institutional_torture_dum[data$education == 3], na.rm = T),
  mean(data$institutional_torture_dum[data$education == 4], na.rm = T),
  mean(data$institutional_torture_dum[data$indigenous == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$indigenous == 0], na.rm = T),
  mean(data$institutional_torture_dum[data$cat_age == 1], na.rm = T),
  mean(data$institutional_torture_dum[data$cat_age == 2], na.rm = T),
  mean(data$institutional_torture_dum[data$cat_age == 3], na.rm = T),
  mean(data$institutional_torture_dum[data$cat_age == 4], na.rm = T),
  mean(data$institutional_torture_dum[data$cat_age == 5], na.rm = T),
  mean(data$institutional_torture_dum[data$cat_age == 6], na.rm = T),
  mean(data$institutional_torture_dum, na.rm = T),
  sum(data$institutional_torture_dum, na.rm = T)
)

col3 <- c(
  mean(data$threat_dum[data$P3_2 == 1], na.rm = T),
  mean(data$threat_dum[data$P3_2 == 2], na.rm = T),
  mean(data$threat_dum[data$P3_2 == 4], na.rm = T),
  mean(data$threat_dum[data$P3_2 == 3], na.rm = T),
  mean(data$threat_dum[data$P3_2 == 5|data$P3_2 == 6], na.rm = T),
  mean(data$threat_dum[data$federal_prison == 0], na.rm = T),
  mean(data$threat_dum[data$federal_prison == 1], na.rm = T),
  mean(data$threat_dum[data$homicide == 1], na.rm = T),
  mean(data$threat_dum[data$rape == 1], na.rm = T),
  mean(data$threat_dum[data$theft == 1], na.rm = T),
  mean(data$threat_dum[data$drug_commerce == 1], na.rm = T),
  mean(data$threat_dum[data$drug_possess == 1], na.rm = T),
  mean(data$threat_dum[data$extortion == 1], na.rm = T),
  mean(data$threat_dum[data$illegal_weapons == 1], na.rm = T),
  mean(data$threat_dum[data$kidnap == 1], na.rm = T),
  mean(data$threat_dum[data$SEXO == 1], na.rm = T),
  mean(data$threat_dum[data$SEXO == 0], na.rm = T),
  mean(data$threat_dum[data$education == 1], na.rm = T),
  mean(data$threat_dum[data$education == 2], na.rm = T),
  mean(data$threat_dum[data$education == 3], na.rm = T),
  mean(data$threat_dum[data$education == 4], na.rm = T),
  mean(data$threat_dum[data$indigenous == 1], na.rm = T),
  mean(data$threat_dum[data$indigenous == 0], na.rm = T),
  mean(data$threat_dum[data$cat_age == 1], na.rm = T),
  mean(data$threat_dum[data$cat_age == 2], na.rm = T),
  mean(data$threat_dum[data$cat_age == 3], na.rm = T),
  mean(data$threat_dum[data$cat_age == 4], na.rm = T),
  mean(data$threat_dum[data$cat_age == 5], na.rm = T),
  mean(data$threat_dum[data$cat_age == 6], na.rm = T),
  mean(data$threat_dum, na.rm = T),
  sum(data$threat_dum, na.rm = T)
)

col4 <- c(
  mean(data$P3_2 == 1, na.rm = T),
  mean(data$P3_2 == 2, na.rm = T),
  mean(data$P3_2 == 4, na.rm = T),
  mean(data$P3_2 == 3, na.rm = T),
  mean(data$P3_2 == 5|data$P3_2 == 6, na.rm = T),
  mean(data$federal_prison == 0, na.rm = T),
  mean(data$federal_prison == 1, na.rm = T),
  mean(data$homicide == 1, na.rm = T),
  mean(data$rape == 1, na.rm = T),
  mean(data$theft == 1, na.rm = T),
  mean(data$drug_commerce == 1, na.rm = T),
  mean(data$drug_possess == 1, na.rm = T),
  mean(data$extortion == 1, na.rm = T),
  mean(data$illegal_weapons == 1, na.rm = T),
  mean(data$kidnap == 1, na.rm = T),
  mean(data$SEXO == 1, na.rm = T),
  mean(data$SEXO == 0, na.rm = T),
  mean(data$education == 1, na.rm = T),
  mean(data$education == 2, na.rm = T),
  mean(data$education == 3, na.rm = T),
  mean(data$education == 4, na.rm = T),
  mean(data$indigenous == 1, na.rm = T),
  mean(data$indigenous == 0, na.rm = T),
  mean(data$cat_age == 1, na.rm = T),
  mean(data$cat_age == 2, na.rm = T),
  mean(data$cat_age == 3, na.rm = T),
  mean(data$cat_age == 4, na.rm = T),
  mean(data$cat_age == 5, na.rm = T),
  mean(data$cat_age == 6, na.rm = T),
  NA,
  NA
)
col5 <- c(
  sum(data$P3_2 == 1, na.rm = T),
  sum(data$P3_2 == 2, na.rm = T),
  sum(data$P3_2 == 4, na.rm = T),
  sum(data$P3_2 == 3, na.rm = T),
  sum(data$P3_2 == 5|data$P3_2 == 6, na.rm = T),
  sum(data$federal_prison == 0, na.rm = T),
  sum(data$federal_prison == 1, na.rm = T),
  sum(data$homicide == 1, na.rm = T),
  sum(data$rape == 1, na.rm = T),
  sum(data$theft == 1, na.rm = T),
  sum(data$drug_commerce == 1, na.rm = T),
  sum(data$drug_possess == 1, na.rm = T),
  sum(data$extortion == 1, na.rm = T),
  sum(data$illegal_weapons == 1, na.rm = T),
  sum(data$kidnap == 1, na.rm = T),
  sum(data$SEXO == 1, na.rm = T),
  sum(data$SEXO == 0, na.rm = T),
  sum(data$education == 1, na.rm = T),
  sum(data$education == 2, na.rm = T),
  sum(data$education == 3, na.rm = T),
  sum(data$education == 4, na.rm = T),
  sum(data$indigenous == 1, na.rm = T),
  sum(data$indigenous == 0, na.rm = T),
  sum(data$cat_age == 1, na.rm = T),
  sum(data$cat_age == 2, na.rm = T),
  sum(data$cat_age == 3, na.rm = T),
  sum(data$cat_age == 4, na.rm = T),
  sum(data$cat_age == 5, na.rm = T),
  sum(data$cat_age == 6, na.rm = T),
  NA,
  NA
) 

table <- cbind(col1, col2, col3, col4, col5)
colnames(table) <- c("Brute Force", "Torture", "Threats", "Proportion", "Total")
rownames(table) <- c(
  "Municipal police",
  "State police",
  "Ministerial police",
  "Federal police",
  "Military",
  "State Prison",
  "Federal Prison",
  "Homicide",
  "Rape",
  "Theft",
  "Drugs (commerce)",
  "Drugs (possession)",
  "Extortion",
  "Illegal weapons",
  "Kidnapping",
  "Male",
  "Female",
  "Primary or less",
  "Secondary",
  "High School",
  "College or postgraduate",
  "Indigenous lang",
  "No indigenous lang",
  "18-25",
  "26-35",
  "36-45",
  "45-55",
  "56-65",
  "65+",
  "Proportion",
  "Number"
)

xtable(
  table, digits = 3
)

###Table: Crimes
# Crimes table ------------------------------------------------------------
#P5_8_1
state_data <- data[data$federal_prison == 0,] %>% data.frame
N <- nrow(state_data)
df <- data.frame(matrix(ncol = 2, nrow = 40))
colnames(df) <- c("full", "subset")
index <- 1
N <- nrow(data)
for(i in 7:25){
  name <- paste0("P5_8_", i)
  df$full[index] <- paste0(round(mean(state_data[[name]], na.rm = TRUE), digits = 4))
  df$full[index + 1] <- paste0("(", 
                               round(
                                 sqrt(1/N*(mean(state_data[[name]], na.rm = T)*(1 - mean(state_data[[name]], na.rm = T)))/N), digits = 4
                               ),
                               ")")
  rownames(df)[c(index, index + 1)] <- c(name, paste(name, "standard error"))
  index <- index + 2
}
name <- "theft"
df$full[39] <- paste0(round(mean(state_data[[name]], na.rm = TRUE), digits = 4))
df$full[40] <- paste0("(", 
                      round(
                        sqrt(1/N*(mean(state_data[[name]], na.rm = T)*(1 - mean(state_data[[name]], na.rm = T)))/N), digits = 6
                      ),
                      ")")

reform <- state_data[!is.na(state_data$allreform),]
index <- 1

N <- nrow(reform)
for(i in 7:25){
  name <- paste0("P5_8_", i)
  df$subset[index] <- paste0(round(mean(reform[[name]], na.rm = TRUE), digits = 4))
  df$subset[index + 1] <- paste0("(", 
                                 round(
                                   sqrt(1/N*(mean(reform[[name]], na.rm = T)*(1 - mean(reform[[name]], na.rm = T)))/N), digits = 4
                                 ),
                                 ")")
  index <- index + 2
}
name <- "theft"
df$subset[39] <- paste0(round(mean(reform[[name]], na.rm = TRUE), digits = 4))
df$subset[40] <- paste0("(", 
                        round(
                          sqrt(1/N*(mean(reform[[name]], na.rm = T)*(1 - mean(reform[[name]], na.rm = T)))/N), digits = 4
                        ),
                        ")")
df %>% xtable()


# Table A5: municipality-month averages  -----------------------------------------------------------------
###aggregate data

vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

datamodels <- 
  data %>% 
  filter(
    !is.na(year_month)&
      !is.na(clave)&
      !is.na(allreform_with_fed)&
      !is.na(group_id)
  ) %>% 
  group_by(
    clave,
    group_id,
    Ent_arrest,
    year_month,
    allreform_with_fed
  ) %>% 
  summarise(
    torture = mean(institutional_torture_dum, na.rm = T),
    force = mean(brute_force_dum, na.rm = T),
    threat = mean(threat_dum, na.rm = T)
  )

x <- lm(torture ~ allreform_with_fed + as.factor(Ent_arrest) + as.factor(year_month),
        datamodels)
x1 <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$torture)]))

x <- lm(force ~ allreform_with_fed + as.factor(Ent_arrest) + as.factor(year_month),
        datamodels)
x2  <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$force)]))

x <- lm(threat ~ allreform_with_fed + as.factor(Ent_arrest) + as.factor(year_month),
        datamodels)
x3 <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$threat)]))


x <- lm(torture ~ allreform_with_fed + as.factor(group_id) + as.factor(year_month),
        datamodels)
x4 <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$torture)]))

x <- lm(force ~ allreform_with_fed + as.factor(group_id) + as.factor(year_month),
        datamodels)
x5  <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$force)]))

x <- lm(threat ~ allreform_with_fed + as.factor(group_id) + as.factor(year_month),
        datamodels)
x6 <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$threat)]))


x <- lm(torture ~ allreform_with_fed + as.factor(clave) + as.factor(year_month),
        datamodels)
x7 <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$torture)]))

x <- lm(force ~ allreform_with_fed + as.factor(clave) + as.factor(year_month),
        datamodels)
x8  <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$force)]))

x <- lm(threat ~ allreform_with_fed + as.factor(clave) + as.factor(year_month),
        datamodels)
x9 <- coeftest(x, vcov = vcovCluster(x, datamodels$group_id[!is.na(datamodels$threat)]))

stargazer(list(x1, x2, x3, x4, x5, x6, x7, x8, x9), 
          omit = c("factor"), star.cutoffs = c(0.05, 0.01, 0.001),
          keep.stat = c('n'))




# Figures A3-5, trends joint operations -------------------------------------------------

###Institutional torture
data %>% 
  filter(
    !is.na(Ent_arrest)
  ) %>% 
  filter(year_arrest > 2000) %>% 
  mutate(
    joint_op_edo = if_else(
      Ent_arrest == 2|Ent_arrest == 8|Ent_arrest == 10|Ent_arrest == 12|Ent_arrest == 16|Ent_arrest == 19|Ent_arrest == 25|Ent_arrest == 28|Ent_arrest == 17|Ent_arrest == 30|Ent_arrest == 15, 1, 0
    )
  ) %>% 
  group_by(
    Ent_arrest
    , joint_op_edo
    , year_arrest
  ) %>% 
  ggplot(
    aes(x = year_arrest, y = institutional_torture_dum)
  ) + stat_smooth(aes(linetype=factor(joint_op_edo), colour=factor(joint_op_edo))) + 
  ggtitle("Trends in Institutionalized Torture around Joint Operations") +
  xlab("Year") + ylab("Torture") + 
  geom_vline(xintercept  = 2006) + 
  labs(color = 'Operation\n', linetype = 'Operation\n') 




data %>% 
  filter(
    !is.na(Ent_arrest)
  ) %>% 
  filter(year_arrest > 2000) %>% 
  mutate(
    joint_op_edo = if_else(
      Ent_arrest == 2|Ent_arrest == 8|Ent_arrest == 10|Ent_arrest == 12|Ent_arrest == 16|Ent_arrest == 19|Ent_arrest == 25|Ent_arrest == 28|Ent_arrest == 17|Ent_arrest == 30|Ent_arrest == 15, 1, 0
    )
  ) %>% 
  group_by(
    Ent_arrest
    , joint_op_edo
    , year_arrest
  ) %>% 
  ggplot(
    aes(x = year_arrest, y = brute_force_dum, colour = factor(joint_op_edo))
  ) + stat_smooth(aes(linetype=factor(joint_op_edo), colour=factor(joint_op_edo))) + 
  ggtitle("Trends in Brute Force around Joint Operations") +
  xlab("Year") + ylab("Brute Force") + 
  geom_vline(xintercept  = 2006) + 
  labs(color = 'Operation\n', linetype = 'Operation\n') 

data %>% 
  filter(
    !is.na(Ent_arrest)
  ) %>% 
  filter(year_arrest > 2000) %>% 
  mutate(
    joint_op_edo = if_else(
      Ent_arrest == 2|Ent_arrest == 8|Ent_arrest == 10|Ent_arrest == 12|Ent_arrest == 16|Ent_arrest == 19|Ent_arrest == 25|Ent_arrest == 28|Ent_arrest == 17|Ent_arrest == 30|Ent_arrest == 15, 1, 0
    )
  ) %>% 
  group_by(
    Ent_arrest
    , joint_op_edo
    , year_arrest
  ) %>% 
  ggplot(
    aes(x = year_arrest, y = threat_dum, colour = factor(joint_op_edo))
  )  + stat_smooth(aes(linetype=factor(joint_op_edo), colour=factor(joint_op_edo))) + 
  ggtitle("Trends in Threats around Joint Operations") +
  xlab("Year") + ylab("Threats") + 
  geom_vline(xintercept  = 2006) + 
  labs(color = 'Operation\n', linetype = 'Operation\n') 



